This script contains the results for all hypotheses using both the MAR and the MNAR imputed datasets for the sensitivity analysis. At the end multiple regression assumptions are checked across the MAR imputed datasets.


Load Packages and Data

library(tidyverse)
library(mice)
library(effsize)
library(rempsyc)
library(car)

imp_MAR <- readRDS("imp_MAR.RData")
imp_MNAR1 <- readRDS("imp_MNAR1.RData")
imp_MNAR2 <- readRDS("imp_MNAR2.RData")
df_final <- readRDS("df_final.RData")

Primary Analysis (MAR)

Hypothesis 1 - Affect measures predict depression retrospectively

m1 <- with(imp_MAR, 
             lm(BDI_prae ~ scale(mean_neg_T5) + scale(mean_pos_T5) + 
                  scale(rmssd_neg_T5) + scale(rmssd_pos_T5)))
nice_table(summary(pool(m1), conf.int = TRUE))

Term

estimate

std.error

statistic

df

p

2.5 %

97.5 %

95% CI

(Intercept)

22.69

0.88

25.82

114.64

< .001***

20.95

24.43

[20.95, 24.43]

scale(mean_neg_T5)

2.23

1.24

1.80

111.54

.075

-0.23

4.68

[-0.23, 4.68]

scale(mean_pos_T5)

2.36

1.45

1.63

108.43

.107

-0.51

5.23

[-0.51, 5.23]

scale(rmssd_neg_T5)

-0.58

1.11

-0.52

115.67

.603

-2.77

1.62

[-2.77, 1.62]

scale(rmssd_pos_T5)

-0.41

1.29

-0.32

114.14

.749

-2.97

2.14

[-2.97, 2.14]

# R-Squared
mean(sapply(m1$analyses, function(m) summary(m)$r.squared))
## [1] 0.03687233

Hypothesis 2 - Affect measures predict depression prospectively

m2 <- with(imp_MAR, 
           lm(BDI_post ~ scale(mean_neg_T55) + scale(mean_pos_T55) + 
                scale(rmssd_neg_T55) + scale(rmssd_pos_T55)))
nice_table(summary(pool(m2), conf.int = TRUE))

Term

estimate

std.error

statistic

df

p

2.5 %

97.5 %

95% CI

(Intercept)

11.88

1.04

11.45

74.60

< .001***

9.82

13.95

[9.82, 13.95]

scale(mean_neg_T55)

-0.84

1.51

-0.56

92.32

.578

-3.83

2.15

[-3.83, 2.15]

scale(mean_pos_T55)

-1.50

1.84

-0.82

79.88

.417

-5.17

2.17

[-5.17, 2.17]

scale(rmssd_neg_T55)

1.52

1.17

1.30

82.02

.199

-0.81

3.85

[-0.81, 3.85]

scale(rmssd_pos_T55)

-0.60

1.53

-0.39

69.29

.697

-3.64

2.45

[-3.64, 2.45]

# R-Squared
mean(sapply(m2$analyses, function(m) summary(m)$r.squared))
## [1] 0.05049228

Hypothesis 3 - Affect measures change over time

Mean Positive Affect

t.test(df_final$mean_pos_T5, df_final$mean_pos_T55, paired = TRUE)
## 
##  Paired t-test
## 
## data:  df_final$mean_pos_T5 and df_final$mean_pos_T55
## t = -1.3065, df = 128, p-value = 0.1937
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.008052211  0.001647616
## sample estimates:
## mean difference 
##    -0.003202297
cohen.d(df_final$mean_pos_T5, df_final$mean_pos_T55, paired = TRUE)
## 
## Cohen's d
## 
## d estimate: -0.1160853 (negligible)
## 95 percent confidence interval:
##       lower       upper 
## -0.29165124  0.05948068

Mean Negative Affect

t.test(df_final$mean_neg_T5, df_final$mean_neg_T55, paired = TRUE)
## 
##  Paired t-test
## 
## data:  df_final$mean_neg_T5 and df_final$mean_neg_T55
## t = 0.87931, df = 128, p-value = 0.3809
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.003986877  0.010364555
## sample estimates:
## mean difference 
##     0.003188839
cohen.d(df_final$mean_neg_T5, df_final$mean_neg_T55, paired = TRUE)
## 
## Cohen's d
## 
## d estimate: 0.07577303 (negligible)
## 95 percent confidence interval:
##       lower       upper 
## -0.09416962  0.24571568

Variation Positive Affect

t.test(df_final$rmssd_pos_T5, df_final$rmssd_pos_T55, paired = TRUE)
## 
##  Paired t-test
## 
## data:  df_final$rmssd_pos_T5 and df_final$rmssd_pos_T55
## t = -0.39184, df = 128, p-value = 0.6958
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.0002449449  0.0001639675
## sample estimates:
## mean difference 
##   -4.048869e-05
cohen.d(df_final$rmssd_pos_T5, df_final$rmssd_pos_T55, paired = TRUE)
## 
## Cohen's d
## 
## d estimate: -0.04176154 (negligible)
## 95 percent confidence interval:
##      lower      upper 
## -0.2517352  0.1682121

Variation Negative Affect

t.test(df_final$rmssd_neg_T5, df_final$rmssd_neg_T55, paired = TRUE)
## 
##  Paired t-test
## 
## data:  df_final$rmssd_neg_T5 and df_final$rmssd_neg_T55
## t = 0.51839, df = 128, p-value = 0.6051
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.0001444791  0.0002470576
## sample estimates:
## mean difference 
##    5.128924e-05
cohen.d(df_final$rmssd_neg_T5, df_final$rmssd_neg_T55, paired = TRUE)
## 
## Cohen's d
## 
## d estimate: 0.05608721 (negligible)
## 95 percent confidence interval:
##      lower      upper 
## -0.1571456  0.2693200

Hypothesis 4 - Change in affect measures predics change in depression

m4 <- with(imp_MAR, 
           lm(delta_BDI ~ scale(delta_mean_neg) + scale(delta_mean_pos) + 
              scale(delta_rmssd_neg) + scale(delta_rmssd_pos)))
nice_table(summary(pool(m4), conf.int = TRUE))

Term

estimate

std.error

statistic

df

p

2.5 %

97.5 %

95% CI

(Intercept)

-10.81

1.23

-8.79

85.17

< .001***

-13.25

-8.36

[-13.25, -8.36]

scale(delta_mean_neg)

2.70

1.84

1.47

90.78

.145

-0.95

6.35

[-0.95, 6.35]

scale(delta_mean_pos)

1.03

2.01

0.51

90.96

.611

-2.97

5.02

[-2.97, 5.02]

scale(delta_rmssd_neg)

-0.95

1.40

-0.68

88.00

.500

-3.73

1.83

[-3.73, 1.83]

scale(delta_rmssd_pos)

0.80

1.61

0.50

81.41

.622

-2.41

4.01

[-2.41, 4.01]

# R-Squared
mean(sapply(m4$analyses, function(m) summary(m)$r.squared))
## [1] 0.03700676

Sensitivity Analysis

MNAR1 Model

Hypothesis 1 - Affect measures predict depression retrospectively

m1_A <- with(imp_MNAR1, 
             lm(BDI_prae ~ scale(mean_neg_T5) + scale(mean_pos_T5) + 
                  scale(rmssd_neg_T5) + scale(rmssd_pos_T5)))
nice_table(summary(pool(m1_A), conf.int = TRUE))

Term

estimate

std.error

statistic

df

p

2.5 %

97.5 %

95% CI

(Intercept)

23.01

0.89

25.88

114.70

< .001***

21.24

24.77

[21.24, 24.77]

scale(mean_neg_T5)

2.36

1.26

1.87

109.24

.065

-0.15

4.86

[-0.15, 4.86]

scale(mean_pos_T5)

2.37

1.45

1.63

111.14

.106

-0.51

5.25

[-0.51, 5.25]

scale(rmssd_neg_T5)

-0.93

1.13

-0.82

112.80

.412

-3.18

1.31

[-3.18, 1.31]

scale(rmssd_pos_T5)

-0.34

1.30

-0.26

115.45

.796

-2.91

2.24

[-2.91, 2.24]

# R-Squared
mean(sapply(m1_A$analyses, function(m) summary(m)$r.squared))
## [1] 0.04062428

Hypothesis 2 - Affect measures predict depression prospectively

m1_A <- with(imp_MNAR1, 
           lm(BDI_post ~ scale(mean_neg_T55) + scale(mean_pos_T55) + 
                scale(rmssd_neg_T55) + scale(rmssd_pos_T55)))
nice_table(summary(pool(m1_A), conf.int = TRUE))

Term

estimate

std.error

statistic

df

p

2.5 %

97.5 %

95% CI

(Intercept)

14.27

1.15

12.44

62.33

< .001***

11.97

16.56

[11.97, 16.56]

scale(mean_neg_T55)

-0.79

1.64

-0.48

80.34

.632

-4.05

2.47

[-4.05, 2.47]

scale(mean_pos_T55)

-1.34

2.06

-0.65

64.18

.518

-5.46

2.78

[-5.46, 2.78]

scale(rmssd_neg_T55)

1.75

1.24

1.41

77.50

.162

-0.72

4.22

[-0.72, 4.22]

scale(rmssd_pos_T55)

-0.70

1.52

-0.46

79.10

.646

-3.72

2.32

[-3.72, 2.32]

# R-Squared
mean(sapply(m1_A$analyses, function(m) summary(m)$r.squared))
## [1] 0.05272962

Hypothesis 4 - Change in affect measures predics change in depression

m4_A <- with(imp_MNAR1, 
           lm(delta_BDI ~ scale(delta_mean_neg) + scale(delta_mean_pos) + 
              scale(delta_rmssd_neg) + scale(delta_rmssd_pos)))
nice_table(summary(pool(m4_A), conf.int = TRUE))

Term

estimate

std.error

statistic

df

p

2.5 %

97.5 %

95% CI

(Intercept)

-8.74

1.31

-6.68

76.57

< .001***

-11.34

-6.13

[-11.34, -6.13]

scale(delta_mean_neg)

2.36

1.95

1.20

81.71

.232

-1.53

6.24

[-1.53, 6.24]

scale(delta_mean_pos)

0.61

2.09

0.29

87.51

.770

-3.55

4.77

[-3.55, 4.77]

scale(delta_rmssd_neg)

-0.92

1.42

-0.65

91.16

.520

-3.74

1.90

[-3.74, 1.90]

scale(delta_rmssd_pos)

0.83

1.61

0.51

88.95

.609

-2.37

4.03

[-2.37, 4.03]

# R-Squared
mean(sapply(m4_A$analyses, function(m) summary(m)$r.squared))
## [1] 0.03104246

MNAR2 Model

Hypothesis 1 - Affect measures predict depression retrospectively

m1_B <- with(imp_MNAR2, 
             lm(BDI_prae ~ scale(mean_neg_T5) + scale(mean_pos_T5) + 
                  scale(rmssd_neg_T5) + scale(rmssd_pos_T5)))
nice_table(summary(pool(m1_B), conf.int = TRUE))

Term

estimate

std.error

statistic

df

p

2.5 %

97.5 %

95% CI

(Intercept)

22.34

0.88

25.42

115.46

< .001***

20.60

24.08

[20.60, 24.08]

scale(mean_neg_T5)

2.03

1.24

1.64

112.77

.104

-0.42

4.48

[-0.42, 4.48]

scale(mean_pos_T5)

2.46

1.42

1.73

115.80

.086

-0.35

5.26

[-0.35, 5.26]

scale(rmssd_neg_T5)

-0.16

1.12

-0.14

112.58

.889

-2.38

2.07

[-2.38, 2.07]

scale(rmssd_pos_T5)

-0.56

1.28

-0.44

117.99

.663

-3.08

1.97

[-3.08, 1.97]

# R-Squared
mean(sapply(m1_B$analyses, function(m) summary(m)$r.squared))
## [1] 0.03495761

Hypothesis 2 - Affect measures predict depression prospectively

m1_B <- with(imp_MNAR2, 
           lm(BDI_post ~ scale(mean_neg_T55) + scale(mean_pos_T55) + 
                scale(rmssd_neg_T55) + scale(rmssd_pos_T55)))
nice_table(summary(pool(m1_B), conf.int = TRUE))

Term

estimate

std.error

statistic

df

p

2.5 %

97.5 %

95% CI

(Intercept)

10.40

1.09

9.52

59.93

< .001***

8.22

12.59

[8.22, 12.59]

scale(mean_neg_T55)

-0.93

1.48

-0.63

90.59

.530

-3.88

2.01

[-3.88, 2.01]

scale(mean_pos_T55)

-1.72

1.78

-0.97

82.86

.337

-5.26

1.82

[-5.26, 1.82]

scale(rmssd_neg_T55)

1.27

1.13

1.12

85.01

.264

-0.98

3.53

[-0.98, 3.53]

scale(rmssd_pos_T55)

-0.43

1.46

-0.30

74.19

.767

-3.34

2.47

[-3.34, 2.47]

# R-Squared
mean(sapply(m1_B$analyses, function(m) summary(m)$r.squared))
## [1] 0.0447545

Hypothesis 4 - Change in affect measures predics change in depression

m4_B <- with(imp_MNAR2, 
           lm(delta_BDI ~ scale(delta_mean_neg) + scale(delta_mean_pos) + 
              scale(delta_rmssd_neg) + scale(delta_rmssd_pos)))
nice_table(summary(pool(m4_B), conf.int = TRUE))

Term

estimate

std.error

statistic

df

p

2.5 %

97.5 %

95% CI

(Intercept)

-11.93

1.27

-9.41

75.23

< .001***

-14.46

-9.41

[-14.46, -9.41]

scale(delta_mean_neg)

2.80

1.77

1.58

98.25

.118

-0.72

6.31

[-0.72, 6.31]

scale(delta_mean_pos)

1.25

2.03

0.62

86.26

.540

-2.78

5.28

[-2.78, 5.28]

scale(delta_rmssd_neg)

-1.01

1.34

-0.76

97.30

.452

-3.67

1.65

[-3.67, 1.65]

scale(delta_rmssd_pos)

0.88

1.61

0.54

79.45

.588

-2.33

4.08

[-2.33, 4.08]

# R-Squared
mean(sapply(m4_B$analyses, function(m) summary(m)$r.squared))
## [1] 0.03698853

Complete Case Analysis

df_cc <- df_final[complete.cases(df_final),]
summary(lm(BDI_prae ~ scale(mean_neg_T5) + scale(mean_pos_T5) + 
                scale(rmssd_neg_T5) + scale(rmssd_pos_T5), data = df_cc))
## 
## Call:
## lm(formula = BDI_prae ~ scale(mean_neg_T5) + scale(mean_pos_T5) + 
##     scale(rmssd_neg_T5) + scale(rmssd_pos_T5), data = df_cc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.0194  -6.0691  -0.7176   5.8495  29.5712 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          22.5068     1.1629  19.354   <2e-16 ***
## scale(mean_neg_T5)    1.0187     1.5503   0.657    0.513    
## scale(mean_pos_T5)    2.8181     1.7066   1.651    0.103    
## scale(rmssd_neg_T5)   0.4251     1.4293   0.297    0.767    
## scale(rmssd_pos_T5)  -0.7623     1.5747  -0.484    0.630    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.936 on 68 degrees of freedom
## Multiple R-squared:  0.04512,    Adjusted R-squared:  -0.01105 
## F-statistic: 0.8034 on 4 and 68 DF,  p-value: 0.5273
summary(lm(BDI_post ~ scale(mean_neg_T55) + scale(mean_pos_T55) + 
                scale(rmssd_neg_T55) + scale(rmssd_pos_T55), data = df_cc))
## 
## Call:
## lm(formula = BDI_post ~ scale(mean_neg_T55) + scale(mean_pos_T55) + 
##     scale(rmssd_neg_T55) + scale(rmssd_pos_T55), data = df_cc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.871  -7.350  -1.913   6.501  25.056 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           12.0993     1.1470  10.548 5.89e-16 ***
## scale(mean_neg_T55)   -1.1343     1.6832  -0.674   0.5026    
## scale(mean_pos_T55)   -0.9809     1.8535  -0.529   0.5984    
## scale(rmssd_neg_T55)   3.1844     1.3608   2.340   0.0222 *  
## scale(rmssd_pos_T55)  -1.6283     1.5039  -1.083   0.2828    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.8 on 68 degrees of freedom
## Multiple R-squared:  0.09975,    Adjusted R-squared:  0.04679 
## F-statistic: 1.884 on 4 and 68 DF,  p-value: 0.1233
summary(lm(delta_BDI ~ scale(delta_mean_neg) + scale(delta_mean_pos) + 
              scale(delta_rmssd_neg) + scale(delta_rmssd_pos), data = df_cc))
## 
## Call:
## lm(formula = delta_BDI ~ scale(delta_mean_neg) + scale(delta_mean_pos) + 
##     scale(delta_rmssd_neg) + scale(delta_rmssd_pos), data = df_cc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.4422  -7.9622  -0.2727   8.5198  29.3597 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -10.4075     1.3647  -7.626 1.04e-10 ***
## scale(delta_mean_neg)    1.6662     1.9776   0.843    0.402    
## scale(delta_mean_pos)    1.0327     2.2195   0.465    0.643    
## scale(delta_rmssd_neg)   0.5902     1.5464   0.382    0.704    
## scale(delta_rmssd_pos)   1.5504     1.7867   0.868    0.389    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.66 on 68 degrees of freedom
## Multiple R-squared:  0.04569,    Adjusted R-squared:  -0.01045 
## F-statistic: 0.8139 on 4 and 68 DF,  p-value: 0.5207

Check multiple regression assumptions (in the MAR datasets)

Homoscedasticity and normality

Model 1

plot_residuals <- function(model, m) {
  results <- model[["analyses"]][[m]]
  residuals <- residuals(results)
  fitted_values <- fitted(results)
  
  par(mfrow = c(1, 2))
  
  #Residual vs Fitted Plot
  plot(fitted_values, residuals, 
       xlab = "Fitted values", ylab = "Residuals", 
       pch = 20, col = "blue")
  abline(h = 0, col = "red", lwd = 2)
  
  #Q-Q Plot
  qqnorm(residuals, main = "")
  qqline(residuals, col = "red", lwd = 2)
  
   mtext(paste("Imputation", m), outer = TRUE,line = -3)
}

for (i in 1:50) {
  plot_residuals(m1, m = i)
}

#### Model 2

for (i in 1:50) {
  plot_residuals(m2, m = i)
}

#### Model 4

for (i in 1:50) {
  plot_residuals(m4, m = i)
}

Multicollinearity

Model 1

completed_df <- complete(imp_MAR)
vif(lm(BDI_prae ~ scale(mean_neg_T5) + scale(mean_pos_T5) + 
                  scale(rmssd_neg_T5) + scale(rmssd_pos_T5), 
       data = completed_df))
##  scale(mean_neg_T5)  scale(mean_pos_T5) scale(rmssd_neg_T5) scale(rmssd_pos_T5) 
##            1.926390            2.585032            1.589765            2.131584

Model 2

vif(lm(BDI_prae ~ scale(mean_neg_T55) + scale(mean_pos_T55) + 
                  scale(rmssd_neg_T55) + scale(rmssd_pos_T55), 
       data = completed_df))
##  scale(mean_neg_T55)  scale(mean_pos_T55) scale(rmssd_neg_T55) 
##             2.408126             3.276685             1.349086 
## scale(rmssd_pos_T55) 
##             2.043889

Model 3

vif(lm(delta_BDI ~ scale(delta_mean_neg) + scale(delta_mean_pos) + 
                   scale(delta_rmssd_neg) + scale(delta_rmssd_pos), 
       data = completed_df))
##  scale(delta_mean_neg)  scale(delta_mean_pos) scale(delta_rmssd_neg) 
##               2.312139               2.774158               1.310595 
## scale(delta_rmssd_pos) 
##               1.657758